1.1 Scaletta

1. Cosa prevedere: i passaggi fondamentali
2. Definizione del problema
3. Raccolta dei dati
4. Scelta di un modello di previsione

1.2 Previsioni

Spesso basata su ipotesi più o meno verificate, alle volte influenzata da distorsioni cognitive latenti, quasi sempre caratterizzata da informazioni scarse o addirittura insufficienti, una previsione in condizioni di incertezza circa le più disparate situazioni è necessaria.

Una previsione è necessaria in tantissime situazioni, sia che venga poi formulata in modo esplicito o in modo implicito, il secondo tende a prevalere. Per il direttore di una filiale, decidere quanti moduli per la sottoscrizione di servizi di investimento stampare ogni mattina comporta una valutazione del numero atteso giornaliero di sottoscrizioni. Per una banca decidere quante nuove filiali aprire nei prossimi cinque anni richiede una previsione della domanda futura di servizi bancari. Sempre per una banca, la pianificazione dei turni di personale per il call centre richiede previsioni dei volumi delle chiamate. Stimare il livello delle sofferenze previste per il trimestre in corso è ugualmente importante, per il management di un importante gruppo bancario.

1.3 Determinare cosa prevedere

Nelle prime fasi di un progetto di previsione, è necessario prendere decisioni su ciò che deve essere previsto. Una volta determinato quali siano le previsioni necessarie, è necessario trovare o raccogliere i dati su cui si baseranno le previsioni. I dati necessari per la previsione potrebbero già esistere. Al giorno d’oggi, molti dati si rendono disponibili in modo più o meno automatico, nonostante il loro impiego possa essere più o meno facile.

Un’attività di previsione prevede in genere tre passaggi fondamentali:

1. La definizione del problema.
2. La raccolta dei dati.
3. La scelta di un modello di previsione.

1.4 Un esempio: gli NPL in Italia

Torniamo alla fine del 2016 —malgrado i dati in nostro possesso siano aggiornati ad oggi. I non-performing loan sono un tema di attualità importante. E preoccupante al tempo stesso.

library(tidyverse)
data <- read_csv2("data/TRI30265.csv", skip = 1)
# Tidy data.
data <- data %>% filter(!is.na(date)) %>% 
  rename(def_npe = defaulted_npe)
library(lubridate)
data <- data %>% 
    select(date, cured_exp, def_npe) %>% 
    mutate(date = dmy(date)) %>% 
    mutate_at(c("cured_exp", "def_npe"), ~(./1e+06))

data <- data %>% mutate(def_npe = def_npe - cured_exp) %>% 
  arrange(date) %>% 
  mutate(def_npe = round(def_npe, 2))
npe <- ts(data$def_npe,
          start = c(1999, 1), 
          freq = 4)
(npe2 <- window(npe, end = c(2016,4)))

#>       Qtr1  Qtr2  Qtr3  Qtr4
#> 1999  71.9  71.6  69.4  69.1
#> 2000  65.6  65.0  64.6  63.8
#> 2001  59.5  56.7  48.9  48.3
#> 2002  49.0  48.5  48.2  49.0
#> 2003  49.3  50.3  51.0  52.8
#> 2004  55.4  55.8  56.2  57.0
#> 2005  57.3  56.7  57.5  48.4
#> 2006  48.3  48.4  48.5  48.8
#> 2007  50.0  49.8  50.8  48.4
#> 2008  50.3  47.3  45.9  41.9
#> 2009  46.1  51.1  57.9  62.5
#> 2010  67.5  71.5  75.6  81.4
#> 2011  86.7  90.7 106.0 109.8
#> 2012 112.7 117.1 122.9 127.3
#> 2013 134.1 141.0 149.2 155.7
#> 2014 164.1 172.4 178.2 177.2
#> 2015 184.2 188.6 196.0 197.2
#> 2016 195.3 197.7 198.2 200.7

library(forecast)
autoplot(npe2) + 
  labs(title = "Gli NPL in Italia al 31/12/2016", 
       caption = "Fonte: tavola TRI30265, Banca d'Italia, https://infostat.bancaditalia.it", 
       x = "tempo (dati trimestrali)", 
       y = "NPL, mld di euro")

fcast_2017_2018 <- forecast(npe2, h = 8)

knitr::kable(fcast_2017_2018, digits = 2)

autoplot(fcast_2017_2018) + 
  labs(title = "Previsioni su NPL in Italia al 31/12/2016", 
       caption = "Fonte: tavola TRI30265, statistical data base, Banca d'Italia, https://infostat.bancaditalia.it", 
       x = "tempo (dati trimestrali)", 
       y = "NPL, mld di euro")

1.5 Take-away

Una previsione, anche se in condizioni di incertezza, circa le più disparate situazioni si rende spesso necessaria. Ma non è banale: prendere decisioni su ciò che deve essere previsto, raccogliere i dati, scegliere un modello di previsione sono passaggi per nulla scontati.

2.1 Scaletta

1. Graficare i dati
2. Tipologie di serie storiche
3. Trend
4. Ciclo
5. Stagionalità
6. Stazionarietà

2.2 Quali serie storiche

La prima cosa da fare, come in qualsiasi attività di analisi dei dati, è quello di tracciare graficamente i dati. I grafici consentono di visualizzare molte caratteristiche, osservazioni insolite, modifiche nel tempo e relazioni tra i dati. Per applicare un modello di previsione delle serie storiche temporali appropriato è necessario individuare le caratteristiche dei dati e il luogo naturale con cui iniziare è quello di utilizzare i grafici.

Per i dati di serie temporali, il grafico ovvio per iniziare è un grafico temporale. Utilizzando, ad esempio, i dati mensili relativi alle vendite al dettaglio non alimentari rilevato dall’ISTAT (https://www4.istat.it/it/prodotti/banche-dati/serie-storiche) e scaricabili da https://drive.google.com/open?id=1bn6fy9VGAsiLZYt-RCEQKkzkgC5EsURi, le osservazioni sono tracciate rispetto al tempo di osservazione, con rilevazioni consecutive unite da linee rette.

library(tidyverse)
dt_ts <- read_csv2("data/vendite_dettaglio.csv", col_names = "vendite")
library(forecast)
retail <- ts(dt_ts$vendite, 
    start = c(2010, 1),
    end = c(2019, 5),
    frequency = 12,
    names = "vendite"
    )

Un grafico rivela immediatamente alcune caratteristiche interessanti.

autoplot(retail) +
  ggtitle("Indice ISTAT delle vendite al dettaglio (non alimentari)") +
  xlab("tempo (dati mensili)") +
  ylab("100 = base indice dicembre 2009")

2.3 Tipologie di serie storiche

Una serie storica macroeconomica può essere descritta in termini di:

• trend (andamento tendenziale); una serie storica ha un trend quando c’è un aumento o una diminuzione di medio/lungo termine nei dati. Non deve essere necessariamente lineare e a volte è chiaramente evidente un suo cambiamento.

autoplot(retail) + 
  stat_smooth(method = "loess", se = FALSE, span = 0.80) + 
  labs(title = "Trend", 
       subtitle = "Indice ISTAT delle vendite al dettaglio (non alimentari)", 
       x = "tempo (dati mensili)", 
       y = "100 = base indice dicembre 2009")

• stagionalità; quando una serie temporale è influenzata da fattori ricorrenti stagionali legati al periodo dell’anno (mese o trimestre) o il giorno della settimana si definisce stagionale.

ggseasonplot(retail, year.labels=TRUE, year.labels.left=TRUE) + 
  labs(title = "Stagionalità", 
       subtitle = "Grafico stagionale (linea per anno)\nIndice ISTAT delle vendite al dettaglio (non alimentari)", 
       x = "tempo (dati mensili)", 
       y = "100 = base indice dicembre 2009")

ggseasonplot(retail, polar = TRUE) + 
  labs(title = "Stagionalità", 
       subtitle = "Grafico stagionale in coord. polari (linea per anno)\nIndice ISTAT vendite al dettaglio (non aliment.)") + 
  scale_color_discrete(name = "anno") + 
  xlab("") + 
    theme(legend.position = "bottom")

• ciclicità; una serie storica è ciclica, e non stagionale, quando l’aumento o la diminuzione dei dati non avviene con frequenza fissa —come in quelle stagionali. Queste fluttuazioni sono di solito dovute a condizioni economiche e sono spesso correlate al “ciclo economico”. La durata di queste fluttuazioni è di solito di almeno 2 anni.

library(fpp2)
autoplot(uschange[, c("Consumption", "Income")]) +
  xlab("tempo (rilevazioni trimestrali)") + ylab("variazione %") + 
  stat_smooth(method = "loess", se = FALSE, span = 0.15) + 
  labs(title = "Ciclicità", 
       subtitle = "Ciclo economico di consumi e reddito negli USA") +
  guides(colour=guide_legend(title = "serie storica")) + 
    theme(legend.position = "bottom")

• stazionarietà; una serie storica è stazionaria quando i suoi valori non dipendono dal momento in cui sono rilevati. A questo modo, una serie storica con un evidente trend (crescente o decrescente) o con una elevata componente stagionale tende ad essere non stazionaria.

2.4 Take-away

L’analisi dei dati comincia con un grafico, il cui scopo è quello di mettere in evidenza degli schemi ricorrenti nelle rilevazioni storiche. Sono soprattutto trend, ciclo e stagionalità che caratterizzano una serie storica, oltre al carattere di stazionarietà.

3.1 Scaletta

1. L’indice di correlazione e la funzione di autocorrelazione (ACF)
2. La funzione di autocorrelazione parziale (PACF)
3. Il White noise

3.2 La funzione di autocorrelazione

Nello stesso modo in cui l’indice di correlazione (\(\rho\)) misura l’entità della relazione (lineare) esistente tra due variabili (reddito e consumi, ad esempio), l’indice di autocorrelazione è una misura della relazione (lineare) tra gli effetti ritardati che si susseguono nella serie temporale. Esistono tanti coefficienti di autocorrelazione quanti “ritardi” (lag) è possibile immaginare

\[\rho_{y_t, y_{t - 1}} \quad \rho_{y_t, y_{t - 2}} \ldots \quad \rho_{y_t, y_{t - n}}\]

Il grafico che segue contiene la serie storica dei dati sulle variazioni percentuali trimestrali dei consumi negli USA dal 1970 al 2016.

library(tidyverse)
library(fpp2)
consum <- uschange[ , "Consumption"]
autoplot(consum) + labs(
  x = "tempo (trimestri)", 
  y = "var. %", 
  title = "Consumi USA")

Ad esempio, si può immaginare che esista una correlazione tra il dato sulla variazione dei consumi trimestrali negli USA di un trimestre rispetto al trimestre precedente (1 ritardo trimestrale); oppure si può ipotizzare, ad esempio, che esista una relazione stagionale per cui una correlazione tra il dato dei consumi di un trimestre e quello dello stesso trimestre nell’anno precedente (4 ritardi trimestrali, ossia un anno). Ma osservando il grafico precedente non è facile notare se esistano o meno questi effetti ritardati tra i dati successivi nei consumi USA.

Per questo si ricorre quasi sempre ad un grafico, denominato ACF per autocorrelation function, che contiene gli indici di correlazione misurati su ritardi successivi crescenti come quello riportato di seguito e relativo ai dati sui consumi trimestrali negli USA.

acf_cons <- ggAcf(consum) + 
  labs(title = "consumi Usa", 
       x = "ritardo (lag)", 
       y = "funzione di autocorr.\n(ACF)")
acf_cons

Il primo segmento verticale nel grafico rappresenta l’indice di correlazione tra le variazioni trimestrali nei consumi USA rispetto al trimestre precedente; il secondo segmento rappresenta l’indice di correlazione tra le variazioni trimestrali nei consumi USA rispetto a 2 trimestri precedenti; e così via. In questo caso, la correlazione tra trimestri successivi (1 ritardo) è circa 0,35, la correlazione tra i dati sui consumi in un semestre (2 trimestri, o 2 ritardi, corrispondo ad un semestre) è circa 0,31, e la correlazione tra i dati sui consumi in un anno (4 trimestri, o 4 ritardi, ossia 1 anno) è circa 0,14.

I coefficienti di autocorrelazione aiutano a identificare dei modelli, degli schemi ricorrenti (pattern) tipici delle serie storiche di tipo stagionale. Nei grafici seguenti sono riportati i dati mensili relativi all’indice delle vendite al dettaglio non alimentari pubblicato dall’Istat e i relativi indici di autocorrelazione (ACF). È facile osservare il “picco” (spike) positivo di correlazione a 12 mesi (quasi 0,9) e a 24 mesi (0,84 circa). Inoltre, l’indice di correlazione a 4 mesi di ritardo è negativo, anche se solo di poco (–0,32 circa).

Le vendite al dettaglio tendono ad essere caratterizzate da una forte componente stagionale; il dato di ogni mese:

1. è fortemente e positivamente correlato con quello dello stesso mese dell’anno precedente;
2. è negativamente correlato con quello di quattro mesi prima.

In presenza di una “forte” autocorrelazione è possibile quindi letteralmente “sfruttare” i dati passati per elaborare delle previsioni sufficientemente accurate circa i dati futuri.

3.3 Funzione di autocorrelazione parziale

Se i dati di una serie storica (\(y_t\)) sono autocorrelati tra di loro, immaginiamo vi sia una autocorrelazione con un ritardo di 1 periodo (1 lag), questo significa che \(y_{t}\) e \(y_{t - 1}\) sono correlati così come \(y_{t-1}\) e \(y_{t - 2}\) sono correlati; questo significa anche che \(y_{t}\) e \(y_{t - 2}\) tenderanno ad essere correlati perché “collegati” entrambi a \(y_{t - 1}\): il grafico di partial autocorrelation function (PACF) riporta i coefficienti di autocorrelazione parziale a ritardi successivi \(t-1\), \(t - 2\), \(t - 3\), \(\ldots\) calcolati eliminando l’effetto delle autocorrelazioni intermedie. Il grafico che segue rappresenta la PACF delle vendite al dettaglio ISTAT.

ggPacf(retail) +
  labs(title = "Forte autocorrelazione a 12 mesi per \nle vendite al dettaglio", 
       x = "ritardo (lag mensile)", 
       y = "funzione di autocorr. parziale\n(PACF)")

La correlazione tra i dati a 24 mesi è condizionata al fatto che i \(y_t\) e \(y_{t-24}\) sono entrambi legati tra di loro da \(y_{t-12}\): questo è evidente nell’osservare nel grafico PACF che l’unico picco rimanente è a 12 mesi, quello a 24 mesi è scomparso; l’autocorrelazione segnala una stagionalità significativa a 12 mesi, non a 24 mesi.

3.4 White noise

Serie storiche che non mostrano autocorrelazioni vengono denominate white noise. Nel grafico che segue un esempio illustre, e largamente atteso: l’S&P 500 (o meglio il suo ETF, i dati sono scaricabili da https://drive.google.com/file/d/1v94zkqNKmFxrGO5BPw2XNI_qVdPpHJUn/view?usp=sharing).

etfs <- read_csv("data/ETFs.csv")
spy <- ts(
  etfs$SPY, 
  start = c(2007, 1, 1), 
  end = c(2019, 6, 30), 
  frequency = 12, 
  names = "SPY"
)

autoplot(spy) + labs(title = "White noise apparente per l'S&P 500",
                  x = "tempo (mesi)", 
                  y = "rendimenti mensili")

Per le serie storiche white noise ci aspettiamo che ogni coefficiente di autocorrelazione sia approssimativamente pari a zero. Chiaramente, non tutti saranno esattamente pari a zero. Per una serie storica white noise ci aspettiamo che il 95% dei coefficienti di correlazione sia all’interno di un intervallo di ampiezza pari a \(\pm1{,}96/\sqrt{T}\), dove \(T\) rappresenta la lunghezza della serie temporale.

Nel caso dell’S&P 500 (o meglio del suo ETF) ci aspettiamo che il 95% degli indici di autocorrelazione cada all’interno di una fascia di ampiezza \(\pm 1{,}96/\sqrt{150}\) ossia circa \(\pm\) 0,16. Nell’autocorrelogramma riportato sopra tutti gli indici di correlazione —ad eccezione di quello a 6 mesi— sono inferiori a \(|\) 0,16 \(|\).

Nei grafici seguenti sono riportate le funzioni di autocorrelazione (ACF) sia dello S&P 500 che dei dati relativi ai consumi trimestrali negli USA.

ggAcf(spy) + 
  labs(title = "S&P 500", 
       x = "ritardo (lag)", 
       y = "funzione di autocorr.\n(ACF)") + 
  coord_cartesian(ylim = c(-0.25, 0.35))

acf_cons + 
  coord_cartesian(ylim = c(-0.25, 0.35))

Si nota chiaramente come la funzione di autocorrelazione (ACF) per l’S&P 500 non sia significativamente diversa da zero (nessun indice, se non quello a 6 mesi) esce dalle bande blu tratteggiate, mentre quella relativa ai consumi trimestrali negli USA è significativamente positiva almeno fino a 4 trimestri di ritardo.

Se da un lato i dati passati relativi ai consumi USA forniscono informazioni utili circa quelli successivi (e quindi futuri), così non è per i dati relativi all’S&P 500. Possiamo quindi concludere che i rendimenti mensili dell’S&P 500 non siano altro che white noise.

Cosa significa che una serie storica (si tratta in particolare di quelle finanziarie) è di tipo white noise? La risposta è esplicitamente contenuta nella nota avvertenza “i rendimenti passati non sono indicativi di quelli futuri”.

3.5 Take-away

L’autocorrelazione in una serie storica è uno dei motori principali di previsione. L’autocorrelazione misura l’entita della relazione esistente tra rilevazioni successive di rilevazioni in una serie storica. In presenza di autocorrelazione la capacità previsiva di un fenomeno è fortemente avvantaggiata. In sua assenza, i dati passati non sono indicativi di quelli futuri.

4.1 Scaletta

1. Descrizione del modello
2. Scelta dell’ordine di autoregressione
3. Esame del grafico ACF e PACF
4. Previsioni con il modello

4.2 Il modello

L’impiego dei modelli autoregressivi ipotizza, ragionevolmente, che la variabile di una serie temporale —il PIL, ad esempio— sia in relazione con i suoi valori passati. Infatti, è abbastanza intuitivo che il passato immediato di una variabile possa possedere una capacità previsiva sufficientemente utile per elaborare delle proiezioni a breve, medio o lungo termine.

Schematizzando con più formalità, si ipotizza che la serie storica di un fenomeno \(y_t\) sia governata da un modello autoregressivo con la seguente formulazione

\[y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \ldots + \phi_p y_{t-p} + \varepsilon_t\]

dove \(\varepsilon_t\) è la componente white noise. Si rappresenta quindi un modello autoregressivo \(y_t\) come funzione di una costante \(c\) e dei \(p\) valori passati della serie storica \(y_{t-1}\), \(y_{t-2}\), \(\ldots\), \(y_{t-p}\)—oltre alla componente residuale di disturbo \(\varepsilon\).

Un fenomeno economico come l’andamento del PIL dell’area euro (riportato nel grafico seguente e i cui dati sono scaricabili da https://drive.google.com/file/d/1CUPXeBI7fgtTJ0Q4_Kqtymk4oSpcAZvQ/view?usp=sharing) può essere rappresentato ai fini previsivi utilizzando un modello autoregressivo i cui coefficienti \(c\), \(\phi_1\), \(\phi_2\), \(\ldots\), \(\phi_p\) devono essere stimati utilizzando i dati passati.

library(tidyverse)
eu_gdp <- read_csv("data/eu_gdp.csv")
eugdp_ts <- ts(eu_gdp$NAEXKP01EZQ659S, 
                start = 1996, 
                frequency = 4, 
                names = "eugdp")

library(forecast)
autoplot(eugdp_ts) + 
  labs(title = "L'andamento del PIL nell'area euro dal 1996", 
       x = "tempo (dati trimestrali)", 
       y = "PIL reale (var. %)"
       )

Se ipotizziamo che il modello autoregressivo sia composto da una sola variabile, ossia sia presente un solo ritardo allora

\[y_t = c + \phi_1 y _{t-1} + \varepsilon_t\]

definiamo il modello come modello autoregressivo di primo ordine, indicandolo con AR(1).

Utilizzando ad esempio i dati sulle variazioni percentuali trimestrali del PIL dell’area euro dal gennaio 1996 riportati nel grafico precedente proviamo a stimare un modello autoregressivo di primo ordine (AR1).

(fit_ar1 <- Arima(eugdp_ts, order = c(1, 0, 0)))

#> Series: eugdp_ts 
#> ARIMA(1,0,0) with non-zero mean 
#> 
#> Coefficients:
#>         ar1   mean
#>       0.899  1.524
#> s.e.  0.042  0.706
#> 
#> sigma^2 estimated as 0.578:  log likelihood=-106
#> AIC=219   AICc=219   BIC=226

Dall’output dell’elaboratore il modello autoregressivo che rappresenta l’andamento temporale del PIL nell’area euro è dato da una costante \(c = 1{,}524 \left(1 - 0{,}899 \right) = 0{,}154\), da un coefficiente \(\phi_1 = 0{,}8987\) e da una componente di disturbo casuale \(\varepsilon_t\) con una deviazione standard pari a \(\sqrt{0{,}578} = 0{,}76\), ossia

\[\begin{split}y_t &= 1{,}524 \times \left(1 - 0{,}899 \right) + 0{,}899 y_{t-1} + \varepsilon_t \\ &= 0{,}154 + 0{,}899 y_{t-1} + \varepsilon_t\end{split}\]

Come si è visto la costante \(c\) non è inclusa nell’output dell’elaboratore ma deve essere ricavata tramite un accordimento \(c = \left(1 - \phi_1 - \phi_2 - \ldots - \phi_p\right)\mu\), dove \(\mu\) è riportato nell’output come mean: questo perché il tipo di parametrizzazione scelto in alcuni ambienti (R in questo caso) è diverso da quello usato nell’equazione generale riportata all’inizio.

La capacità del modello, che abbiamo appena stimato, di rappresentare il fenomeno è indicata dalle statistiche AIC, AICc e BIC (quanto minori, tanto migliore il modello); scegliendo la seconda (AICc, Akaike Information Coefficient corrected) proviamo ora confrontare il modello AR(1) con un secondo modello di secondo ordine AR(2) la cui stima è rappresentata di seguito.

(fit_ar2 <- Arima(eugdp_ts, order = c(2, 0, 0)))

#> Series: eugdp_ts 
#> ARIMA(2,0,0) with non-zero mean 
#> 
#> Coefficients:
#>         ar1     ar2   mean
#>       1.501  -0.657  1.564
#> s.e.  0.075   0.075  0.364
#> 
#> sigma^2 estimated as 0.322:  log likelihood=-79.1
#> AIC=166   AICc=167   BIC=176

In questo caso, il modello AR(2) può essere riformulato nel modo seguente

\[\begin{split} y_t &= \left(1 - 1{,}501 + 0{,}657\right) \times 1{,}564 + 1{,}501 y_{t-1} - 0{,}657 y_{t-2} + \varepsilon_t \\ &= -0{,}24 + 1{,}501 y_{t-1} - 0{,}657 y_{t-2} + \varepsilon_t \end{split}\]

Di fatto questo modello AR(2) possiede una capacità di adattamento ai dati (e quindi anche previsiva) maggiore di quello AR(1): la statistica AICc è pari a 167 che rispetto a quella del modello AR(1), pari a 219, è sensibilmente inferiore. Proviamo quindi ad utilizzare questo secondo modello AR(2) per fare delle previsioni sul PIL dell’area euro a due anni, ossia 8 trimestri.

fit_ar2 %>% forecast(h = 8) %>% 
  autoplot() + 
  labs(title = "previsioni su 2 anni (8 trimestri) si un modello AR(2)", 
       x = "tempo (dati trimestrali)", 
       y = "PIL reale (var. %)")

Le bande azzurro scuro e azzurro chiaro sono gli intervalli di previsione rispettivamente all’80% e al 95% di probabilità e rappresentano l’entità dell’incertezza relativa alla previsione basata sul modello che è stato stimato. I modelli autoregressivi sono molto utilizzati nel campo delle previsioni economiche e soprattutto macroeconomiche perché si prestano ad una facile integrazione con i modelli moving average nei diffusissimi modelli ARIMA.

Normalmente, la scelta dei numero di ritardi (o variabili) in modello autoregressivo è affidata all’esame del grafico della funzione di autocorrelazione, il cui campione nei dati rilevati può fornire degli elementi utili per la scelta. Per un modello autoregressivo \(AR(p)\) la scelta di \(p\) è effettuata a questo modo: se i grafici ACF e PACF dei dati mostrano questi pattern:

• il grafico ACF decade esponenzialmente o in modo sinusoidale;
• c’è un picco significativo nel grafico PACF al ritardo \(p\) e successivamente più nessuno.

ggAcf(eugdp_ts) + 
  labs(title = "Decadimento sinusoidale di ACF")
ggPacf(eugdp_ts) + 
  labs(title = "2 spike nei primi lag di PACF")

In questo caso, ad esempio, il fatto che vi siano 4 spike significativi nei primi 4 ritardi nella funzione ACF e solo 2 spike significativi nella funzione PACF suggerisce un modello AR(2), come quello selezionato da noi.

Chiaramente, la scelta del modello può essere affidata —come normalmente è— ad algoritmi con maggiore capacità computazionale rispetto all’esame visivo dei grafici di ACF e PACF.

4.3 Take-away

Nei modelli autoregressivi la variabile di una serie temporale —–il PIL, ad esempio–— è in relazione con i suoi valori passati; e questa relazione viene sfrutta ai fini previsivi. La scelta dell’ordine di autoregressione — il numero di lag— può essere facilitata dall’esame visivo dei grafici di ACF e PACF.

5.1 Scaletta

1. Descrizione del modello
2. Indici di bontà di adattamento (AIC, AICc, BIC)
3. Previsioni con il modello

5.2 Il modello

Un modello moving average utilizza gli errori nelle previsioni passate (shock) come variabili previsive in un modo simile a quello di un modello di regressione multipla, con l’importante differenza che le variabili rappresentate dagli errori non sono effettivamente rilevate, ma devono essere stimate con delle procedure iterative.

Un modello moving average di \(y_t\) è funzione di una costante \(c\), dei \(q\) errori nelle previsioni passate e di una componente residuale di disturbo \(\varepsilon_t\), di seguito

\[y_t = c + \varepsilon_t + \theta_1 \varepsilon_{t-1} + \theta_2 \varepsilon_{t-2} + \ldots + \theta_q \varepsilon_{t-q}\]

dove \(\varepsilon_t\) è la componente di disturbo white noise. In un modello moving average la serie temporale \(y_t\) è rappresentata da una media mobile degli errori passati. Di per sé, un modello moving average non è niente altro che un modello autoregressivo infinito, è solo che per motivi di parsimonia computazionale si preferisce rappresentarlo diversamente.

Un modello moving average non deve essere confuso con lo smoothing moving average: il primo è impiegato per fare delle previsioni, mentre il secondo per stimare la componente trend-ciclo dei valori passati.

L’andamento del PIL dell’area euro (riportato nel grafico seguente, i cui dati sono scaricabili da https://drive.google.com/file/d/1CUPXeBI7fgtTJ0Q4_Kqtymk4oSpcAZvQ/view?usp=sharing) può essere rappresentato ai fini previsivi utilizzando un modello moving average i cui coefficienti \(c\), \(\theta_1\), \(\theta_2\), \(\ldots\), \(\theta_q\) devono essere stimati utilizzando i dati a disposizione.

library(tidyverse)
eu_gdp <- read_csv("data/eu_gdp.csv")
eugdp_ts <- ts(eu_gdp$NAEXKP01EZQ659S, 
                start = 1996, 
                frequency = 4, 
                names = "eugdp")

library(forecast)
autoplot(eugdp_ts) + 
  labs(title = "L'andamento del PIL nell'area euro dal 1996", 
       x = "tempo (dati trimestrali)", 
       y = "PIL reale (var. %)"
       )

Ipotizzando che il modello moving average sia composto da una sola variabile, ossia

\[y_t = c + \varepsilon_t + \theta_1 \varepsilon_{t-1}\]

proviamo a stimare i parametri \(\theta_1\) e \(c\) del modello moving average che abbiamo ipotizzato con i dati relativi alle variazioni trimestrali del PIL nell’area euro dal 1996.

(fit_ma1 <- Arima(eugdp_ts, order = c(0, 0, 1)))

#> Series: eugdp_ts 
#> ARIMA(0,0,1) with non-zero mean 
#> 
#> Coefficients:
#>         ma1   mean
#>       1.000  1.578
#> s.e.  0.038  0.196
#> 
#> sigma^2 estimated as 0.922:  log likelihood=-129
#> AIC=265   AICc=265   BIC=272

Il modello che abbiamo appena stimato può quindi essere rappresentato così

\[y_t = 1{,}578 + \varepsilon_t + 1{,}000 \varepsilon_{t-1}\]

L’algoritmo, più o meno standard, per la stima di un modello moving average —ma anche per i modelli AR, così come per i più generali modelli ARIMA— è il metodo di massima verosimiglianza (Maximum Likelihood Estimation, o MLE) in cui i parametri restituiti —in questo caso \(\theta_1\), \(\theta_2\), \(\ldots\)— corrispondono a quelli che rendono massima la probabilità di ottenere, con il modello ipotizzato, i dati effettivamente osservati.

Nell’output riportato sopra viene restituito il logaritmo della probabilità di ottenere i dati osservati con i parametri del modello ipotizzato (log likelihood): evitando inutili dettagli, il miglior modello è quello che massimiizza questo dato, in altre parole è il modello che meglio si adatta ai dati osservati.

In pratica, la capacità del modello di adattarsi ai dati e quindi di essere utile per fare delle previsioni può essere rappresentata tramite le statistiche AIC, AICc e BIC (quanto minori, tanto migliore il modello), che consentono di confrontare la bontà descrittiva e previsiva di due o più modelli statistici.

Scegliendo la seconda (AICc, Akaike Information Coefficient corrected) proviamo ora confrontare il modello MA(1) con un secondo modello di secondo ordine MA(2) la cui stima è rappresentata di seguito.

(fit_ma2 <- Arima(eugdp_ts, order = c(0, 0, 2)))

#> Series: eugdp_ts 
#> ARIMA(0,0,2) with non-zero mean 
#> 
#> Coefficients:
#>         ma1    ma2   mean
#>       1.466  0.507  1.573
#> s.e.  0.073  0.066  0.234
#> 
#> sigma^2 estimated as 0.606:  log likelihood=-109
#> AIC=226   AICc=226   BIC=236

Questo secondo modello MA(2) può formulato allora in questo modo

\[y_t = 1{,}573 + \varepsilon_t + 1{,}466 \varepsilon_{t-1} + 0{,}507 \varepsilon_{t-2}\]

La statistica AICc del secondo modello MA(2) indica che il primo modello MA(1) è meno capace di adattarsi ai dati, e per questo motivo dovremmo utilizzare il modello MA(2) per fare delle previsioni circa l’evoluzione futura del PIL nell’area euro.

fit_ma2 %>% forecast(h = 8) %>% 
  autoplot() + 
  labs(title = "previsioni su 2 anni (8 trimestri) di un modello MA(2)", 
       x = "tempo (dati trimestrali)", 
       y = "PIL reale (var. %)")

Le bande azzurro scuro e azzurro chiaro sono gli intervalli di previsione rispettivamente all’80% e al 95% di probabilità e rappresentano l’entità dell’incertezza relativa alla previsione basata sul modello che è stato stimato. I modelli moving average sono molto utilizzati nel campo delle previsioni economiche e soprattutto macroeconomiche perché si prestano ad una facile integrazione con i modelli autoregressivi nei diffusissimi modelli ARIMA.

5.3 Take-away

Un modello moving average utilizza gli shock passati come variabili previsive in un modo simile a quello di un modello di regressione multipla, con l’importante differenza che le variabili rappresentate dagli errori (o shock) non sono effettivamente rilevate, ma devono essere stimate con delle procedure iterative. La bontà di un modello —e in questo caso di un modello moving average— è giudicata in funzione di alcuni indici (AIC, AICc o BIC).

6.1 Scaletta

1. Condizione di stazionarietà
2. Differencing
3. Unit root test (URT)

6.2 La stazionarietà in una serie storica

Una serie storica è stazionaria se le sue proprietà non dipendono dal momento temporale in cui viene osservata. Per questo motivo le serie temporali con evidenti trend o con caratteri di stagionalità altrettanto evidenti non possono essere considerate stazionarie. Una serie temporale stazionaria non si caratterizza per schemi ricorrenti prevedibili nel tempo —ancora, trend o stagionalità. Una serie storica stazionaria ha un grafico temporale orizzontale che si caratterizza per una variabilità costante.

library(tidyverse)

eu_gdp <- read_csv("data/eu_gdp.csv")
eugdp_ts <- ts(eu_gdp$NAEXKP01EZQ659S, 
                start = 1996, 
                frequency = 4, 
                names = "eugdp")

it_retail <- read_csv2("data/vendite_dettaglio.csv", col_names = "vendite")
it_retail_ts <- ts(it_retail$vendite, 
    start = c(2010, 1),
    end = c(2019, 5),
    frequency = 12,
    names = "vendite")

Ad esempio, l’andamento del PIL nell’area euro dal 1996, riportato nel grafico di sinistra seguente (i dati sono scaricabili da questo link https://drive.google.com/file/d/1CUPXeBI7fgtTJ0Q4_Kqtymk4oSpcAZvQ/view?usp=sharing) è una serie storica stazionaria, mentre l’indice delle vendite al dettaglio (non alimentari) ISTAT riportato nel grafico di destra (i dati sono disponibili tramite questo link https://drive.google.com/open?id=1bn6fy9VGAsiLZYt-RCEQKkzkgC5EsURi) non è una serie storica stazionaria per l’evidente tratto di stagionalità e per l’apparente trend decrescente riscontrabile nel periodo tra il 2010 e il 2014.

library(forecast)

autoplot(eugdp_ts) + 
  labs(title = "L'andamento del PIL \nnell'area euro dal 1996", 
       x = "tempo (dati trimestrali)", 
       y = "PIL reale (var. %)")

autoplot(it_retail_ts) +
  ggtitle("Indice ISTAT delle \nvendite al dettaglio (non alimentari)") +
  xlab("tempo (dati mensili)") +
  ylab("100 = base indice dicembre 2009")

Perché ci interessa che una serie temporale sia stazionaria? Perché è indispensabile per applicare un modello previsivo. Se una serie storica non è stazionaria si può intervenire per renderla stazionaria tramite una trasformazione, di per sé, molto semplice: calcolando le differenze tra osservazioni successive. Questa trasformazione è chiamata differencing.

library(tidyquant)
aapl_prices  <- tq_get("AAPL", 
                       get = "stock.prices", 
                       from = " 2010-01-01") %>% 
  tq_transmute(select = adjusted, 
               mutate_fun = to.monthly)

aapl_prices_ts <- ts(aapl_prices$adjusted, 
                     start = c(2010, 1), 
                     frequency = 12, 
                     names = "AAPL")

In tante situazioni il differencing consente di trasformare una serie storica non stazionaria in una che lo sia e su cui si possa applicare un modello previsivo. Ad esempio, la serie storica mensile del prezzo di chiusura di Apple dal gennaio 2010, riportata nel grafico seguente sulla sinistra si caratterizza per un evidente trend crescente.

autoplot(aapl_prices_ts) +
  ggtitle("Prezzo di chiusura di \nApple dal 2010") +
  xlab("tempo (dati mensili)") +
  ylab("AAPL \n(adjusted close)")

autoplot(diff(aapl_prices_ts, differences = 1)) +
  ggtitle("Differencing nel prezzo di \nchiusura di Apple dal 2010") +
  xlab("tempo (dati mensili)") +
  ylab("differencing di AAPL \n(adjusted close)")

Nel grafico di destra è riportata la serie storica dopo il differencing (che in questo caso non è altro che di fatto la variazione di prezzo del titolo Apple), si nota immediatamente che il trend crescente è stato eliminato. Ma la serie storica sembra non ancora stazionaria a causa della variabilità delle variazioni di prezzo che tende ad aumentare nel tempo (questo è dovuto allo stesso trend crescente presente nei dati iniziali).

6.3 Unit root test

Il primo esame per verificare se una serie temporale è stazionaria è quello visivo, mediante l’osservazione del grafico dei dati. Un secondo esame più formale, e quindi più rigoroso, che consente di determinare se la serie storica debba essere sottoposta al differencing, perché non stazionaria, è il test di radice unitaria (unit root test).

library(urca)
ur.kpss(eugdp_ts) %>% summary()

#> 
#> ####################### 
#> # KPSS Unit Root Test # 
#> ####################### 
#> 
#> Test is of type: mu with 3 lags. 
#> 
#> Value of test-statistic is: 0.323 
#> 
#> Critical value for a significance level of: 
#>                 10pct  5pct 2.5pct  1pct
#> critical values 0.347 0.463  0.574 0.739

Lo unit root test di Kwiatkowski-Phillips-Schmidt-Shin (KPSS, riportato sopra) applicato alle serie storica dei dati sul PIL nell’area euro restituisce una statistica pari a 0,323: affinchè la serie storica sia stazionaria è necessario che la statistica dello unit root test di KPSS sia inferiore ai valori critici (0,347, 0,463, …). In questo caso quindi, si può concludere che la serie storica del PIL dell’area euro sia stazionaria.

ur.kpss(aapl_prices_ts) %>% summary()

#> 
#> ####################### 
#> # KPSS Unit Root Test # 
#> ####################### 
#> 
#> Test is of type: mu with 4 lags. 
#> 
#> Value of test-statistic is: 2.2 
#> 
#> Critical value for a significance level of: 
#>                 10pct  5pct 2.5pct  1pct
#> critical values 0.347 0.463  0.574 0.739

Lo stesso test KPSS applicato alla serie storica dei dati di Apple (di sopra) invece restituisce un dato molto maggiore (2,2) dei valori critici e per questo motivo non si può considerare la serie come stazionaria. Se sottraiamo i valori successivi (differencing) della serie otteniamo invece

ur.kpss(diff(aapl_prices_ts)) %>% summary()

#> 
#> ####################### 
#> # KPSS Unit Root Test # 
#> ####################### 
#> 
#> Test is of type: mu with 4 lags. 
#> 
#> Value of test-statistic is: 0.0773 
#> 
#> Critical value for a significance level of: 
#>                 10pct  5pct 2.5pct  1pct
#> critical values 0.347 0.463  0.574 0.739

La serie storica di Apple dopo il differencing è invece stazionaria, applicando il test KPSS la statistica è pari a 0,0866, inferiore al valore critico.

6.4 Take-away

Perché ci interessa che una serie temporale sia stazionaria? Perché è indispensabile per applicare un modello previsivo. Se una serie storica non è stazionaria si può intervenire per renderla stazionaria tramite una trasformazione, di per sé, molto semplice: calcolando le differenze tra osservazioni successive. Questa trasformazione è chiamata differencing.

7.1 Scaletta

1. Descrizione del modello ARIMA
2. I grafici ACF e PACF
3. Gli algoritmi automatici

7.2 Il modello

Cosa succede se integriamo i modelli autoregressivi (AR) con quelli moving average (MA)? Otteniamo i modelli ARMA. Per poter applicare un modello previsionale è necessario che la serie storica sia stazionaria, se non lo è si ricorre al differencing per tentare di indurre la stazionarietà. E se integriamo i modelli AR con quelli MA aver operato le dovute trasformazioni di differencing per rendere stazionaria la serie temporale? Otteniamo una famiglia di modelli, non stagionali, denominata AutoRegressive Integrated Moving Average, ma meglio conosciuta con l’acronimo ARIMA. Con la notazione consueta un modello ARIMA è rappresentato così

\[y'_{t} = c + \phi_1y'_{t-1} + \ldots + \phi_p y'_{t-p} + \theta_1 \varepsilon_{t-1} + \ldots + \theta_q \varepsilon_{t - q}+\varepsilon_t\]

dove con \(y'_{t}\) si intende la serie storica dopo il differencing operato \(d\) volte, con \(\phi_p\) il parametro di autoregressione di ordine \(p\), con \(\theta_q\) il parametro di moving average di ordine \(q\). Un modello ARIMA viene descritto utilizzando la notazione \(ARIMA\left(p, d, q\right)\).

Ad esempio, proviamo a considerare la serie storica delle variazioni trimestrali del PIL nell’area euro (i dati sono scaricabili da questo link https://drive.google.com/file/d/1CUPXeBI7fgtTJ0Q4_Kqtymk4oSpcAZvQ/view?usp=sharing) riportata di seguito.

library(tidyverse)
eu_gdp <- read_csv("data/eu_gdp.csv")
eugdp_ts <- ts(eu_gdp$NAEXKP01EZQ659S, 
                start = 1996, 
                frequency = 4, 
                names = "eugdp")

library(forecast)
autoplot(eugdp_ts) + 
  labs(title = "L'andamento del PIL \nnell'area euro dal 1996", 
       x = "tempo (dati trimestrali)", 
       y = "PIL reale (var. %)")

7.3 Esame dei grafici ACF e PACF

La scelta dell’ordine degli ordini \(p\) e \(q\), alle volte, può essere effettuata tramite l’esame visivo dei grafici ACF e PACF nel caso in cui si ipotizzino modelli ARIMA del tipo \(\left(p, d, 0\right)\) o \(\left(0, d, q\right)\), assumendo di fatto che si tratti di un modello autoregressivo AR e di un modello moving average MA nel secondo caso.

I dati provengono da un modello \(\left(p, d, 0\right)\) se i grafici ACF e PACF dei dati (già differenziati) mostrano questi modelli:

• il grafico ACF decade esponenzialmente o in modo sinusoidale;
• c’è un picco significativo nel grafico PACF al ritardo \(p\) e successivamente più nessuno.

I dati provengono da un modello \(\left(0, d, q\right)\) se i grafici ACF e PACF dei dati (già differenziati) mostrano questi modelli:

• il grafico PACF decade esponenzialmente o in modo sinusoidale;
• c’è un picco significativo nel grafico ACF al ritardo \(q\) e successivamente più nessuno.

ggAcf(eugdp_ts) +
  labs(title = "", 
       x = "ritardo (lag trimestrale)", 
       y = "funzione di autocorr. \n(ACF)")
ggPacf(eugdp_ts) +
  labs(title = "", 
       x = "ritardo (lag trimestrale)", 
       y = "funzione di autocorr. parziale\n(PACF)")

Dall’esame dei grafici ACF e PACF per il PIL nell’area euro si può ipotizzare un modello \(ARIMA\left(2, 0, 0\right)\): i picchi nel grafico ACF sono a 1, 2, 3 e 4 lag (ritardi) (probabilmente), ed è evidente una tendenza a decadere in modo sinusoidale, nel grafico PACF notiamo i picchi ai lag 1 e 2.

\[\text{PIL}_t = c + \phi_1\text{PIL}_{t-1} + \phi_2\text{PIL}_{t-2} +\varepsilon_t\]

(fit_arima202 <- Arima(eugdp_ts, order = c(2, 0, 0)))

#> Series: eugdp_ts 
#> ARIMA(2,0,0) with non-zero mean 
#> 
#> Coefficients:
#>         ar1     ar2   mean
#>       1.501  -0.657  1.564
#> s.e.  0.075   0.075  0.364
#> 
#> sigma^2 estimated as 0.322:  log likelihood=-79.1
#> AIC=166   AICc=167   BIC=176

Il modello \(ARIMA\left(2, 0, 0\right)\) stimato può essere riformulato nel modo seguente

\[\begin{split}\text{PIL}_t &= \left(1 - 1{,}501 + 0{,}657 \right) \times 1{,}564 + 1{,}501 \cdot \text{PIL}_{t-1} + \\ & - 0{,}657 \cdot \text{PIL}_{t-2} +\varepsilon_t = \\ & = 0{,}244 + 1{,}501 \cdot \text{PIL}_{t-1} - 0{,}657 \cdot \text{PIL}_{t-2} +\varepsilon_t \end{split} \]

La statistica di bontà nell’adattamento ai dati (goodness-of-fit) AICc è pari a 167.

7.4 La scelta del modello ARIMA e gli algoritmi automatici

Proviamo ora a vedere un modello ARIMA generato da un algoritmo automatico.

(fit_autoarima <- auto.arima(eugdp_ts, 
                             seasonal = FALSE))

#> Series: eugdp_ts 
#> ARIMA(5,0,1) with non-zero mean 
#> 
#> Coefficients:
#>         ar1    ar2    ar3     ar4    ar5    ma1   mean
#>       0.688  0.599  -0.41  -0.368  0.241  0.941  1.568
#> s.e.  0.106  0.129   0.13   0.122  0.102  0.056  0.397
#> 
#> sigma^2 estimated as 0.273:  log likelihood=-70.3
#> AIC=157   AICc=158   BIC=177

Si tratta di un modello \(ARIMA\left(5, 0, 1\right)\) la cui statistica AICc è pari a 158, leggermente migliore del modello \(ARIMA\left(2, 0, 0\right)\) ipotizzato da noi inizialmente e sulla base dei grafici ACF e PACF. Utilizziamo questo modello per fare delle previsioni sul PIL nell’area euro a due anni (8 trimestri).

forecast(fit_autoarima, h = 8) %>% 
  autoplot() + labs(
    title = "Previsioni con ARIMA(5, 0, 1) a 8 trimestri (2 anni)", 
    x = "tempo (trimestri)", 
    y = "var. trimestrali (%) PIL area euro")

La procedura “manuale” corretta per la scelta del modello ARIMA può essere descritta a questo modo:

1. Graficare i dati ed esaminare se vi siano delle osservazioni inconsuete;
2. Se necessario trasformare i dati (Box-Cox, ad esempio) per stabilizzare la variabilità;
3. Se i dati non sono stazionari, effettuare dei differencing (quasi sempre 1 ordine di differencing è sufficiente);
4. Esaminare i grafici ACF/PACF: si tratta di un \(ARIMA(p, d, 0)\) o di un \(ARIMA(0, d, q)\)?
5. Testare il modello prescelto confrontandolo con altri tramite l’AICc per individuarne uno migliore
6. Controllare i residui per verificare che si tratti di white noise (no autocorrelazione);
7. Una volta che i residui sono white noise, generare le preferenze.

Attualmente, è diffusa la prassi —fortemente raccomandata— di impiegare degli algoritmi automatici che consentano la selezione del corretto ordine \(p\), \(d\), e \(q\) dei modelli ARIMA. L’algoritmo automatico di Hyndman-Khandakar (del pacchetto forecast in R) procede nel modo seguente:

1. Il numero di differenze \(d\) (con \(0\leq d \leq2\)) viene determinato utilizzando test KPSS ripetuti.
2. I valori di \(p\) e \(q\) vengono quindi scelti minimizzando l’AICc dopo aver differenziato i dati \(d\) volte. Anziché considerare ogni possibile combinazione di \(p\) e \(q\), l’algoritmo utilizza una ricerca graduale per attraversare l’insieme dei modelli disponibili.
   1. Sono disponibili quattro modelli iniziali:
      1. ARIMA(0,0,0),
      2. ARIMA(2,2.2),
      3. ARIMA(1,0,0),
      4. ARIMA(0,d.1). Una costante è inclusa a meno che non sia \(d = 2\). Se \(d \leq 1\), viene “adattato” un modello aggiuntivo: ARIMA(0,d,0) senza costante.
   2. Il miglior modello (con il valore AICc più piccolo) adattato nel passo (a) è impostato per essere il “modello attuale”.
   3. Si considerano le variazioni del modello attuale:
      1. variare \(p\) e/o \(q\) dal modello attuale di \(\pm 1\);
      2. includere/escludere \(c\) dal modello corrente. Il miglior modello considerato finora (il modello corrente o una di queste varianti) diventa il nuovo modello attuale.
   4. Ripetere il passaggio precedente fino a quando non viene trovato l’AICc inferiore.

7.5 Take-away

I modelli ARIMA sono una famiglia di modelli, non stagionali, che consente con ampia flessibilità di adattare un modello quantitativo ai dati passati di una serie storica per poter effettuare previsioni accurate. Dall’esame dei grafici ACF e PACF è possibile guidare la scelta degli ordini autoregressivi e di moving average di un modello ARIMA.

8.1 Scaletta

1. Scomposizione
2. Exponential smoothing
3. Metodo di Holt-Winters

8.2 La scomposizione

Una serie storica di dati macroeconomici si caratterizza spesso, ma non sempre naturalmente, per l’evidente presenza di tre componenti: (1) trend, (2) ciclo e (3) stagionalità.

library(tidyverse)
dt_ts <- read_csv2("data/vendite_dettaglio.csv", col_names = "vendite")
library(forecast)
retail <- ts(dt_ts$vendite, 
    start = c(2010, 1),
    end = c(2019, 5),
    frequency = 12,
    names = "vendite"
    )

Se prendiamo in esame, nel grafico che segue ad esempio, le rilevazioni mensili dell’indice relativo alle vendite al dettaglio (non alimentari) dell’ISTAT (disponibile da questo link https://drive.google.com/open?id=1bn6fy9VGAsiLZYt-RCEQKkzkgC5EsURi) è facile osservare l’ultima componente, quella stagionale: con una ricorrenza annua fissa e su base mensile si notano evidenti movimenti congiunti nell’indice.

autoplot(retail) +
  ggtitle("Indice ISTAT delle vendite al dettaglio (non alimentari)") +
  xlab("tempo (dati mensili)") +
  ylab("100 = base indice dicembre 2009")

È chiaro anche come certe variazioni dell’indice siano ragionevolmente prevedibili: al forte picco verso l’alto di dicembre segue una contrazione prima a gennaio e poi a febbraio, ad esempio. Questa forte componente stagionale può nascondere ad esempio le altre due, quella di tendenza (trend) e quella ciclica. Infatti, può sembrare che nel periodo che va tra il 2010 e il 2014 fosse presente un trend decrescente nell’indice ISTAT delle vendite al dettaglio. Ma non è molto evidente.

La scomposizione è “organizzata”, per così dire, aggregando le componenti trend e ciclo in un’unica componente: sono le variazioni ricorrenti del trend, ma non ad intervalli regolari su base annua, che descrivono l’andamento ciclico.

8.3 I metodi di scomposizione

Esistono diversi metodi di scomposizione di una serie storica come quello classico (che forma le basi di scomposizione degli altri), l’X11, il SEATS e l’STL. Da un punto di vista formale nel metodo classico, operiamo normalmente una decomposizione della serie storica secondo due approcci, denominati rispettivamente additivo e moltiplicativo. Secondo il primo le rilevazioni della serie storica sono date dalla somma di:

\[y_t = T_t + S_t + R_t\]

dove con \(y_t\), \(T_t\), \(S_t\) e \(R_t\) si indicano la serie storica, la componente trend-ciclo, la componente stagionale e la componente residuale. Il metodo moltiplicativo —che fornisce delle utili proprietà matematiche— è dato da

\[y_t = T_t \times S_t \times R_t\]

Se applichiamo il metodo di scomposizione classico di tipo moltiplicativo possiamo rappresentare l’indice delle vendite al dettaglio ISTAT come è riportato nel grafico seguente.

retail %>% 
  decompose(type = "multiplicative") %>% 
  autoplot() + 
  labs(title = "Scomposizione classica (moltipl.) dell'indice ISTAT", 
       x = "tempo (rilevazioni mensili)")

A questo punto, oltre all’evidente componente stagionale è stata per così dire “filtrata” la componente di trend decrescente che, effettivamente, ha caratterizzato l’andamento delle vendite al dettaglio tra il 2010 e il 2014. Uno dei difetti del metodo di scomposizione classico è, ad esempio, l’assenza della componente di trend nei primi e ultimi 6 mesi (nel caso di serie storiche mensili).

8.4 L’exponential smoothing

Nei modelli di exponential smoothing, un approccio alternativo di rappresentazione di una serie storica caratterizzata da trend e andamento stagionale, è quello di ricorrere ad una forma per componenti in cui la serie storica viene distinta in 4 elementi (sempre distinguendo tra metodo additivo e moltiplicativo: (1) la previsione, (2) la componente di livello (o di smoothing), (3) la componente di trend e (4) la componente stagionale. Il metodo più conosciuto per catturare trend e stagionalità nei modelli di exponential smoothing è quello di Holt-Winters che applicato alla serie storica (sia nella “versione” additiva sia in quella moltiplicativa) produce la previsione per l’indice delle vendite al dettaglio ISTAT riportata nel grafico seguente.

fit1 <- hw(retail, seasonal = "additive")
fit2 <- hw(retail, seasonal = "multiplicative")
autoplot(retail) +
  autolayer(fit1, series="HW previsione (additivo)", PI = FALSE) +
  autolayer(fit2, series = "HW previsione (moltiplicativo)",
    PI = FALSE) +
  xlab("Year") +
  ylab("indice ISTAT vendite dettaglio") +
  ggtitle("Previsione con metodo Holt-Winters") +
  guides(colour=guide_legend(title = "Previsione")) + 
    theme(legend.position = "bottom")

Il metodo additivo è da preferirsi quando le variazioni stagionali sono approssimativamente costanti, mentre è opportuno optare per quello moltiplicativo quando le variazioni stagionali cambiano in modo proporzionale al livello della serie.

8.5 Take-away

Una serie storica di dati macroeconomici si caratterizza spesso, ma non sempre naturalmente, per l’evidente presenza di tre componenti: (1) trend, (2) ciclo e (3) stagionalità. Quest’ultima può nascondere le altre due. Per questo motivo le serie storiche vengono scomposte con lo scopo di mettere in evidenza certe caratteristiche rispetto ad altre.

9.1 Scaletta

1. Scomposizione in parti
2. Media mobile
3. Serie destagionalizzata
4. Metodi X11, SEATS e STL

9.2 La scomposizione in parti

Le serie temporali si presentano con varietà differenti di schemi ricorrenti (pattern) e spesso può essere utile distinguere quelle “parti”, che rappresentino ciascuna uno schema, tramite una scomposizione in sottoelementi, o componenti. Se consideriamo nel grafico che segue l’indice delle vendite al dettaglio (non alimentari) ISTAT, rilevato con frequenza mensile e disponibile da questo link https://drive.google.com/open?id=1bn6fy9VGAsiLZYt-RCEQKkzkgC5EsURi, è molto evidente la componente stagionale.

library(tidyverse)
dt_ts <- read_csv2("data/vendite_dettaglio.csv", col_names = "vendite")

retail <- ts(dt_ts$vendite, 
    start = c(2010, 1),
    end = c(2019, 5),
    frequency = 12,
    names = "vendite"
    )

library(forecast)

autoplot(retail) +
  ggtitle("Indice ISTAT delle vendite al dettaglio (non alimentari)") +
  xlab("tempo (dati mensili)") +
  ylab("100 = base indice dicembre 2009")

9.3 La serie storica destagionalizzata

Se rimuoviamo la componente stagionale dai dati originali rimaniamo con una serie storica destagionalizzata. Anche se la scomposizione di una serie storica viene usata per descrivere i dati storici, può essere anche usata ai fini previsivi. Da un punto di vista didattico il metodo di scomposizione più utilizzato è quello classico. Questo metodo di scomposizione può essere di tipo additivo o di tipo moltiplicativo e, ipotizzando una scomposizione moltiplicativa, la serie storica può essere rappresentata così

\[y_t = \hat{S}_t \times \hat{T}_t \times \hat{R}_t\]

dove \(\hat{S}_t\), \(\hat{T}_t\) e \(\hat{R}_t\) sono rispettivamente stime delle componenti stagionale, di trend e residuale.

Rimuovendo la componente stagionale \(S_t\) dalla serie storica, come si è accennato, rimane la serie storica destagionalizzata

\[y_t = \hat{S}_t \times \hat{A}_t\]

dove \(\hat{A}_t = \hat{T}_t \times \hat{R}_t\) è la componente destagionalizzata che viene estratta separatamente dopo aver stimato la componente stagionale \(\hat{S}_t\).

Il primo passo per in una scomposizione classica è quello di impiegare una media mobile per stimare la componente di trend (vedi grafico di seguito).

autoplot(retail) +
  autolayer(ma(retail, 5), series = "ma-5") +
  autolayer(ma(retail, 25), series = "ma-25") + 
  autolayer(ma(retail, 37), series = "ma-37") + 
  ggtitle("Media mobile indice ISTAT") +
  xlab("tempo (dati mensili)") +
  ylab("100 = base indice dicembre 2009") 

La media mobile può essere riscritta, con la consueta formulazione, nel modo seguente

\[\text{ma}_m=\frac{1}{m}\sum_{j = -k}^ky_{t+j}\]

dove con \(m = 2k+1\) e \(\text{ma}_m\) sta per media mobile a \(m\) periodi. È evidente che alla stima della componente di trend tramite una media mobile mancano le prime e le ultime \(k\) rilevazioni; di conseguenza non è possibile stimare la componente stagionale per i primi e gli ultimi periodi. Applicando il metodo di scomposizione classico all’indice delle vendite al dettaglio ISTAT, otteniamo il grafico della serie storica destagionalizzata \(y_t/S_t = A_t\) riportato di seguito.

retail %>% 
  decompose(type = "multiplicative") %>% 
  seasadj() %>% 
  autoplot() + 
  labs(
    title = "Indice ISTAT vendite al dettaglio (non alimentari)", 
    subtitle = "Componente destagionalizzata (metodo classico moltiplicativo)", 
    x = "tempo (rilevazioni trimestrali)", 
    y = "indice")

9.4 I metodi di destagionalizzazione

Esistono attualmente diversi metodi di scomposizione e di destagionalizzazione migliori rispetto a quello classico. Tra questi troviamo il metodo X11. X11 è un metodo —sviluppato dal US Census Bureau e da Statistics Canada— per la destagionalizzazione delle serie storiche trimestrali e mensili: basato sulla scomposizione classica abbina però altre caratteristiche che consentono di superare alcuni limiti del metodo “classico” (come ad esempio l’assenza di stime di trend-ciclo per le prime e le ultime rilevazioni).

library(seasonal)
fit_x11 <- retail %>% 
  seas(x11 = "")
autoplot(fit_x11) + labs(title = "Scomposizione X11 dell'indice ISTAT")

La serie storica destagionalizzata dell’indice ISTAT dei dati sulle vendite al dettaglio (non alimentari) è riportata nel grafico seguente assieme alla componente originale e a quella di trend stimate con il metodo X11.

autoplot(retail, series="indice ISTAT") +
  autolayer(trendcycle(fit_x11), series="trend") +
  autolayer(seasadj(fit_x11), series="indice destagionalizzato X11") +
  xlab("tempo (rilevazioni mensili)") + ylab("indice ISTAT vendite al dettaglio") +
  ggtitle("X11: dati destagionalizzati vendite al dettaglio ISTAT") +
  scale_colour_manual(values=c("blue","gray","red"),
             breaks=c("indice ISTAT","indice destagionalizzato X11","trend")) +
    theme(legend.position = "bottom") + 
  guides(colour=guide_legend(title = "serie"))

Un altro metodo per la scomposizione delle serie storiche è SEATS. SEATS (per Seasonal Extraction in ARIMA Time Series) è stato sviluppato dalla Banca di Spagna e ora ampiamente impiegato da diverse agenzie governative a livello mondiale e funziona solamente con rilevazioni trimestrali o mensili.

fit_seats <- retail %>% seas()

La serie storica dell’indice ISTAT destagionalizzata con il metodo SEATS —riportata nel grafico di seguito— è molto simile a quella ottenuta con il metodo X11.

autoplot(retail, series="indice ISTAT") +
  autolayer(trendcycle(fit_seats), series="trend") +
  autolayer(seasadj(fit_seats), series="indice destagionalizzato SEATS") +
  xlab("tempo (rilevazioni mensili)") + ylab("indice ISTAT vendite al dettaglio") +
  ggtitle("SEATS: dati destagionalizzati vendite al dettaglio ISTAT") +
  scale_colour_manual(values=c("blue","gray","red"),
             breaks=c("indice ISTAT","indice destagionalizzato SEATS","trend")) +
    theme(legend.position="bottom") + 
  guides(colour=guide_legend(title = "serie"))

Infine, STL (per Seasonal and Trend decomposition using Loess) è un metodo robusto di scomposizione delle serie storiche con molti vantaggi rispetto ai metodi classico, X11 e SEATS: gestisce qualunque tipo di stagionalità (non solo quelle mensile e trimestrale), la componente stagionale può cambiare nel tempo e può essre controllata dall’utente così come il grado di “levigatura” (o smoothness) della componente trend-ciclica ed è resiliente rispetto alle rilevazioni inusuali (o outlier).

fit_stl <- stl(retail, t.window = 13, 
           s.window = "periodic",
           robust=TRUE)

La scomposizione STL dell’indice ISTAT è riportata di seguito.

fit_stl %>%
  autoplot() + 
  labs(title = "Scomposizione STL dell'indice ISTAT", 
       x = "tempo (rilevazioni mensili)")

Ipotizzando che la componente stagionale rimanga costante (o che cambi molto lentamente nel tempo), per effettuare delle previsioni riguardanti la componente destagionalizzata si applicano i metodi di previsione non stagionale più diffusi. A questo punto, ad esempio, possiamo effettuare una previsione dei dati destagionalizzati come segue.

fit_stl %>% seasadj() %>% naive() %>%
  autoplot() + 
  labs(title = "Indice ISTAT vendite al dettaglio (non alimentari)",
       subtitle = "Previsione naïf della componente destagionalizzata (metodo di scomposizione STL)", 
       x = "tempo (rilevazioni mensili)", 
       y = "indice")

Infine, sempre utilizzando la destagionalizzazione STL sull’indice ISTAT effettuiamo una previsione naïf sui dati destagionalizzati per poi ristagionalizzarli ed effettuare una previsione sulla serie storica originale.

fit_stl %>% forecast(method = "naive") %>% autoplot() + 
  labs(title = "Indice ISTAT vendite al dettaglio (non alimentari)",
       subtitle = "Previsione STL + random walk", 
       x = "tempo (rilevazioni mensili)", 
       y = "indice")

9.5 Take-away

Rimuovendo la componente stagionale dai dati originali rimaniamo con una serie storica destagionalizzata. Per scomporre una serie storica si ricorre agli approcci additivo o moltiplicativo e per destagionalizzarla si usano i metodi X11, SEATS o STL.

10.1 Scaletta

1. Relazioni bi-direzionali simmetriche
2. Descrizione del modello
3. Condizione di stazionarietà
4. Scelta dei lag
5. Previsioni con il modello VAR

10.2 Il modello VAR

Consumption and income are interrelated, don’t you forget about it! ammoniva il prof. Samuelson un suo ex-studente diventato consigliere economico del presidente USA. Da un punto di vista macroeconomico esiste una relazione “bidirezionale” tra consumi (\(y_1\)) e reddito (\(y_2\)): i consumi sono influenzati dal reddito e viceversa. Se siamo interessati a creare delle previsioni riguardanti \(y_1\) partendo dalla serie storica di \(y_2\) tenendo conto dell’effetto che questi hanno su \(y_1\) possiamo fare affidamento sui modelli autoregressivi vettoriali (VAR, per vector autoregressive).

library(forecast)
library(fpp2)
autoplot(uschange[,c("Consumption","Income")]) +
  ylab("var. %") + xlab("tempo (rilevazioni trimestrali)") + 
  ggtitle("Consumi e PIL USA dal 1970 al 2016") + 
  guides(colour=guide_legend(title = "Dati"))

Nel contesto dei modelli VAR tutte le variabili sono trattate in modo simmetrico, vengono cioè “modellate” come se tutte si influenzassero reciprocamente. Utilizzando una notazione specifica al contesto di influenza simmetrica reciproca possiamo descrivere un modello VAR che consenta di effettuare delle previsioni tenendo conto della bidirezionalità a cui abbiamo accennato nel modo seguente

\[\begin{aligned} y_{1,t} &= c_1 + \phi_{11, 1} y_{1, t -1 } + \phi_{12, 1} y_{2, t - 1} + \varepsilon_{1, t} \\ y_{2,t} &= c_2 + \phi_{21, 1} y_{1, t -1 } + \phi_{22, 1} y_{2, t - 1} + \varepsilon_{2, t}\end{aligned}\]

dove \(\varepsilon_{1, t}\) e \(\varepsilon_{2, t}\) sono processi white noise che possono essere contemporaneamente correlati tra loro. Il coefficiente \(\phi_{ii, l}\) cattura l’influenza della variabile \(i\) con se stessa al ritardo \(l\), mentre il coefficiente \(\phi_{ij, l}\) indica l’effetto che la variabile \(i\) ha sulla variabile \(j\) al ritardo \(l\).

Se le serie storiche sono stazionarie, adattiamo dei modelli VAR direttamente alle osservazioni (“VAR in livello”), se invece le serie storiche non sono stazionarie, operiamo un differencing adattando il modello VAR sulle differenze consecutive (“VAR in differenza”). L’adattamento è effettuato usando la tecnica dei minimi quadrati (normalmente).

Esiste una terza possibilità. Che le serie storiche non siano stazionarie ma siano cointegrate: esiste una combinazione lineare che le rende stazionarie. In questo caso viene inclusa una specificazione del modello VAR che include un meccanismo di correzione degli errori (error correction mechanism, a cui ci si riferisce con vector error correction mechanism, o VECM).

10.3 La scelta dei lag

Due scelte devono essere effettuate. Quante variabili (\(y_1\) e \(y_2\), per reddito e consumi, ad esempio) e quanti ritardi (lag). Indicando le prime con \(k\) e i secondi con \(p\), i coefficienti da stimare \(k + pk^2\). Mantenendo basso, in generale, il numero di variabili per evitare che l’errore di stima aumenti troppo utilizziamo degli information criteria per determinare quanti lag debbano essere inclusi nel modello. Riprendendo l’esempio iniziale sui consumi e sul reddito (2 variabili, quindi) otteniamo come output quello che segue.

library(vars)
VARselect(uschange[,c("Consumption","Income")], lag.max=8,
  type="const")[["selection"]]

#> AIC(n)  HQ(n)  SC(n) FPE(n) 
#>      5      1      1      5

• AICc è il corrected Akaike Information Coefficient, in base al quale 5 lag dovrebbero essere utilizzati;
• HQ sta per criterio di Hannan-Quinn: 1 solo lag;
• SC sta per Schwartz Criterion o molto più conosciuto com BIC (per Bayesian Information Coefficient): 1 solo lag;
• FPE, per Final Prediction Error: ancora 5 lag.

Il criterio da preferire nel caso dei modelli VAR (evitando inutili dettagli ovviamente, come ad esempio il fatto che \(AICc\) tende a sovrastimare il numero di lag, ecc…) tende ad essere il BIC o SC, in questo caso.

Optiamo inizialmente per adattare (to fit in inglese) un modello VAR quindi con 1 solo lag, poi proviamo ad incrementare di 1 il numero di lag a salire verso 5 (AICc e FPE indicano 5 lag per il VAR). Come per, ad esempio, i modelli ARIMA, un test Portmanteau di autocorrelazione consente di verificare se vi sia correlazione seriale nei residui (e per questo informazione persa che può essere recuperata aumentando il numero di lag).

var1 <- VAR(uschange[, c("Consumption","Income")], p = 1, type = "const")
serial.test(var1, lags.pt = 10, type = "PT.asymptotic")

#> 
#>  Portmanteau Test (asymptotic)
#> 
#> data:  Residuals of VAR object var1
#> Chi-squared = 49, df = 36, p-value = 0.07

C’è ancora correlazione seriale con 1 solo lag (il p-value è pari a 0,0714, non si può escludere l’ipotesi di autocorrelazione). Aggiungiamo un secondo lag (\(p = 2\)).

var2 <- VAR(uschange[, c("Consumption","Income")], p = 2, type = "const")
serial.test(var2, lags.pt = 10, type = "PT.asymptotic")

#> 
#>  Portmanteau Test (asymptotic)
#> 
#> data:  Residuals of VAR object var2
#> Chi-squared = 48, df = 32, p-value = 0.04

Ancora correlazione seriale (p-value = 0,0363). Aggiungiamo un terzo lag.

var3 <- VAR(uschange[, c("Consumption","Income")], p = 3, type = "const")
serial.test(var3, lags.pt = 10, type = "PT.asymptotic")

#> 
#>  Portmanteau Test (asymptotic)
#> 
#> data:  Residuals of VAR object var3
#> Chi-squared = 34, df = 28, p-value = 0.2

Con \(p = 3\) —ossia con 3 lag— possiamo quindi escludere la possibilità di correlazione seriale tra i residui (p-value = 0,214).

I coefficienti \(\phi_{ij_l}\) e le due costanti per le due equazioni (abbiamo due variabili —\(k = 2\)) sono riportati nell’output R di seguito.

var3

#> 
#> VAR Estimation Results:
#> ======================= 
#> 
#> Estimated coefficients for equation Consumption: 
#> ================================================ 
#> Call:
#> Consumption = Consumption.l1 + Income.l1 + Consumption.l2 + Income.l2 + Consumption.l3 + Income.l3 + const 
#> 
#> Consumption.l1      Income.l1 Consumption.l2      Income.l2 Consumption.l3 
#>         0.1910         0.0784         0.1595        -0.0271         0.2265 
#>      Income.l3          const 
#>        -0.0145         0.2908 
#> 
#> 
#> Estimated coefficients for equation Income: 
#> =========================================== 
#> Call:
#> Income = Consumption.l1 + Income.l1 + Consumption.l2 + Income.l2 + Consumption.l3 + Income.l3 + const 
#> 
#> Consumption.l1      Income.l1 Consumption.l2      Income.l2 Consumption.l3 
#>         0.4535        -0.2730         0.0217        -0.0900         0.3538 
#>      Income.l3          const 
#>        -0.0538         0.3875

10.4 La previsioni con il modello VAR(3)

Le previsioni generate dal modello VAR(3) possono essere rappresentate quindi nel grafico seguente —le fasce in azzurro rappresentano gli intervalli di previsione all’80% (fascia più scura) e al 95% (fascia più chiara).

forecast(var3, h = 8) %>%
  autoplot() + labs(x = "tempo (rilevazioni trimestrali)", 
                    title = "Consumi e reddito negli USA", 
                    subtitle = "Previsioni su 8 trimestri con modello VAR")

10.5 Take-away

Se siamo interessati a creare delle previsioni riguardanti \(y_1\) partendo dalla serie storica di \(y_2\) tenendo conto dell’effetto che questi hanno su \(y_1\) possiamo fare affidamento sui modelli autoregressivi vettoriali (VAR, per vector autoregressive). Per la scelta dei lag in un modello VAR si può ricorrere all’indice BIC. Le previsioni in un modello VAR sono molteplici.